// Perspective projection in row major form
function perspective(fovy, aspect, near, far)
{
	var f = 1.0 / Math.tan(fovy / 2.0);
	
	return [
		f/aspect,		0,			0,			0,
		0, 				f, 			0, 			0,
		0, 				0, 			(far+near)/(near-far), (2.0*far*near)/(near-far),
		0,				0,			-1,			0
		];
}

function mmul(left, right, dim) // square matrix multiplication in row major form
{
	var mat = new Array(dim*dim);
	for(var row = 0; row < dim; ++row)
	{
		for(var col = 0; col < dim; ++col)
		{
			var dot = 0;
			for(var idx = 0; idx < dim; ++idx) dot += left[row * dim + idx] * right[col + idx * dim];
			mat[row * dim + col] = dot;
		}
	}
	
	return mat;
}

function rotx(angle)
{
	var cs = Math.cos(angle), sn = Math.sin(angle);
	return [
		1,		0, 		0,		0,
		0,		cs,		-sn,	0,
		0,		sn,		cs,		0,
		0,		0,		0,		1
		]
			
}

function roty(angle)
{
	var cs = Math.cos(angle), sn = Math.sin(angle);
	return [
		cs,		0, 		sn,		0,
		0,		1,		0,		0,
		-sn,	0,		cs,		0,
		0,		0,		0,		1
		]
			
}

function vdot(a, b, dim)
{
    var sum = 0;
    for(var i = 0; i < dim; ++i) sum += a[i] * b[i];
    return sum;
}

function magn(vec, dim)
{
    return Math.sqrt(vdot(vec, vec, dim));
}

function unit(vec, dim)
{
    var m = magn(vec, dim);
    if(m < 1.0e-7) return vec; // zero mangitude, bad divisor
    return svmul(1.0 / m, vec, dim);
}

function translate(x, y, z)
{
	return [1, 0, 0, x, 0, 1, 0, y, 0, 0, 1, z, 0, 0, 0, 1];
}

function orbit(lon, lat, dist)
{
	return mmul(
		translate(0, 0, -dist),
		mmul(rotx(lat), roty(-lon), 4),
		4);
}

function identity(dim)
{
    var mat = new Array(dim * dim);
    for(var col = 0; col < dim; ++col) 
        for(var row = 0; row < dim; ++row)
            mat[col + row * dim] = (col == row) ? 1.0 : 0.0;
    return mat;
}

function svmul(s, vec, dim)
{
    var result = new Array(dim);
    for(var i = 0; i < dim; ++i) result[i] = s * vec[i];
    return result;
}

function vadd(a, b, dim)
{
    var result = new Array(dim);
    for(var i = 0; i < dim; ++i) result[i] = a[i] + b[i];
    return result;
}

function vsub(a, b, dim)
{
    var result = new Array(dim);
    for(var i = 0; i < dim; ++i) result[i] = a[i] - b[i];
    return result;
}

function mvmul(mat, vec, dim)
{
    var result = new Array(4);
    for(var row = 0; row < dim; ++row)
    {
        var sum = 0;
        for(var col = 0; col < dim; ++col) sum += mat[row * dim + col] * vec[col];
        result[row] = sum; 
    }
    return result;
}

function inverse(input, dim) // gauss-jordan inversion
{
    var left = new Array(dim * dim);
    var right = identity(dim);

    for(var i = 0; i < (dim*dim); ++i) left[i] = input[i];

    function rowadd(src, dest, factor)
    {
        for(var i = 0; i < dim; ++i)
        {
            left[dest * dim + i] += factor * left[src * dim + i];
            right[dest * dim + i] += factor * right[src * dim + i];
        }
    }
    
    for(var row = 0; row < dim-1; ++row) // up->down
    {
        for(var pivot = row+1; pivot < dim; ++pivot)
        {
            if(Math.abs(left[row * dim + row]) < 1.0e-7) continue;
            var factor = -left[pivot * dim + row] / left[row * dim + row];
            rowadd(row, pivot, factor);
        }
    }

    for(var row = dim-1; row > 0; --row) // down->up
    {
        for(var pivot = row-1; pivot >= 0; --pivot)
        {
            if(Math.abs(left[row * dim + row]) < 1.0e-7) continue;        
            var factor = -left[pivot * dim + row] / left[row * dim + row];
            rowadd(row, pivot, factor);
        }
    }
    
    function rowmul(row, factor)
    {
        for(var i = 0; i < dim; ++i) right[row * dim + i] *= factor;
    }

    for(var pivot = 0; pivot < dim; ++pivot) // scaling
    {
        if(Math.abs(left[pivot * dim + pivot]) < 1.0e-7) continue;
        rowmul(pivot, 1.0 / left[pivot * dim + pivot]);
    }

    return right;
}


